function HCP_LRR_predictable_behavior(FC_file, LRR_dir, maxLRR_iter, Nperm, test_metric, intrim_csv, outmat, ...
    restricted_csv, subj_ls, bhvr_ls, colloq_ls)

    % HCP_LRR_predictable_behavior(FC_file, LRR_dir, maxLRR_iter, Nperm, test_metric, intrim_csv, outmat, ...
    %     restricted_csv, subj_ls, bhvr_ls, colloq_ls)
    %
    % Permute behavioral scores across subjects, considering the family structures (multi-level block 
    % permutation). Muilti-level block permutation requires the FSL PALM package.
    % Run elastic net (LRR) on the permuted scores to test the predictability of the original 
    % KRR model.
    %
    % Inputs:
    %   - LRR_dir
    %     The directory which contains the elastic net trained models and testing results.
    %   - maxLRR_iter
    %     Maximal random seed used to split the training-test folds for performing elastic net, e.g. 400.
    %   - Nperm
    %     Number of premutations.
    %   - test_metric
    %     The accuracy metric used to perform statistical testing. Choose from 'predictive_COD' and 'corr'.
    %   - intrim_csv
    %     Full path of the intermediate csv file generated by hcp2block2 function.
    %   - outmat
    %     Full path of the output mat file storing the behaviors whose actual prediction accuracy was significantly
    %     above chance.
    %   - restricted_csv 
    %     Full path of the HCP restricted csv file, containing the family structure information. 
    %   - subj_ls
    %     Full subject list (absolute path).
    %   - bhvr_ls 
    %     List of behavioral measures for which matched AA and matched WA can be found (absolute path).
    %   - colloq_ls
    %     List of colloquial names of behavioral variables. The colloquial names should correspond to the
    %     behavioral names in "bhvr_ls".

    addpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'palm', 'palm-alpha109'))
    addpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'glmnet', 'glmnet_matlab'))

    iter = 40;
    
    %% load shared files
    [bhvr_nm, nbhvr] = CBIG_text2cell(bhvr_ls);
    disp(bhvr_nm)
    disp(nbhvr)
    [colloq_nm, ncolloq] = CBIG_text2cell(colloq_ls);
    if(ncolloq ~= nbhvr)
        error('Number of behavioral names is not equal to number of colloquial names.')
    end
    [subjects, nsub] = CBIG_text2cell(subj_ls);
    alpha_FDR = 0.05;
    
    %% create permutations
    % 1. use `hcp2blocks.m` to generate block definitions
    % 2. use `palm_tree.m` to create exchangeable tree
    % 3. use `palm_permtree.m` to generate permutation indices
    idx_perm_out = fullfile(LRR_dir, 'Pset.mat');
    if(~exist(idx_perm_out, 'file'))
        B = hcp2blocks(restricted_csv, intrim_csv, false, dlmread(subj_ls));
        Ptree = palm_tree(B);
        Pset = palm_permtree(Ptree, Nperm+1, false, true);
        save(idx_perm_out, 'Pset', 'B')
    else
        load(idx_perm_out)
    end
    
    %% load and reshape RSFC 
    load(FC_file)
    features = reshape3d_to_2d(corr_mat);

    %% Read covariates which need to be regressed out from RSFC
    cov_X = load(fullfile(LRR_dir, 'cov_X_58behaviors.mat'));

    %% calculate permuted elastic net accuracies
    metrics = {'corr','COD','predictive_COD','MAE','MAE_norm','MSE','MSE_norm'};
    for b = 1:nbhvr
        LRR_perm(features, cov_X.covariates, b, LRR_dir, maxLRR_iter, Nperm, Pset, bhvr_nm, metrics)
    end

    %% calculate p value for each behavior
    p_perm = zeros(nbhvr, 1);
    avg_stats = [];
    avg_null_stats = []
    for b = 1:nbhvr
        curr_avg_stats = [];
        curr_avg_null_stats = [];
        for i = 1:maxKRR_iter
            opt_fname = fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'results', 'optimal_acc', ...
                [bhvr_nm{b} '_final_acc.mat']);
            if(~exist(opt_fname, 'file'))
                continue
            end
            opt = load(opt_fname);
            Nfolds = length(opt.original_statistics);
            orig_stats = zeros(Nfolds, 1);
            for f = 1:Nfolds
                orig_stats(f) = opt.optimal_statistics{f}.(test_metric);
            end

            acc_out = fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'perm.mat');
            load(acc_out)
    
            curr_avg_stats = cat(1, curr_avg_stats, mean(orig_stats, 1));
            curr_avg_null_stats = cat(1, curr_avg_null_stats, mean(stats_perm.(test_metric), 1));
        end
    
        avg_stats = cat(2, avg_stats, curr_avg_stats);
        avg_null_stats = cat(3, avg_null_stats, curr_avg_null_stats);
        p_perm(b) = length(find( mean(curr_avg_null_stats,1) > mean(curr_avg_stats,1) )) ./ Nperm;
    end

    %% multiple comparisons correction
    H = FDR(p_perm, alpha_FDR);
    sig_perm_idx = sort(H);
    sig_behaviors = bhvr_nm(sig_perm_idx);
    
    switch test_metric
    case 'predictive_COD'
        sig_COD = sig_perm_idx;
        p_COD_perm = p_perm;
        save(outmat, 'p_COD_perm', 'H', 'sig_COD', 'sig_behaviors', 'avg_stats', 'avg_null_stats');
    case 'corr'
        sig_corr = sig_perm_idx;
        p_corr_perm = p_perm;
        save(outmat, 'p_corr_perm', 'H', 'sig_corr', 'sig_behaviors', 'avg_stats', 'avg_null_stats');
    end

    rmpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'palm', 'palm-alpha109'))
    rmpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'glmnet', 'glmnet_matlab'))
    
end

function LRR_perm(features, cov_X, b, LRR_dir, maxLRR_iter, Nperm, Pset, bhvr_nm, metrics)

    % LRR_perm(features, cov_X, b, LRR_dir, maxLRR_iter, Nperm, Pset, bhvr_nm, metrics)
    %
    % 

    for i = 1:maxLRR_iter
        subdir = fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b});
        if(~exist(subdir, 'dir'))
            continue
        end

        fprintf('seed: %d, %s\n', i, bhvr_nm{b})
        load(fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
            ['no_relative_10_fold_sub_list_' bhvr_nm{b} '.mat']));
        Nfolds = length(sub_fold);

        acc_out = fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'perm.mat');
        if(exist(acc_out, 'file'))
            continue
        end

        for f = 1:Nfolds
            test_ind = sub_fold(f).fold_index==1;
            train_ind = ~test_ind;

            %% load cross-validatedly regressed behavioral scores
            y_reg = load(fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'y', ...
                ['fold_' num2str(f)], ['y_regress_' bhvr_nm{b} '.mat']));

            %% split training vs test data for RSFC
            feat_train = features(:, train_ind);
            feat_test = features(:, test_ind);

            %% do confound regression from features, if necessary
            if(~isempty(cov_X) && ~strcmpi(cov_X, 'none'))
                [feat_train, beta] = CBIG_regress_X_from_y_train(feat_train', ...
                    cov_X(train_ind, :));
                beta_pre = load(fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
                    'params', ['fold_' num2str(f)], 'feature_regress_beta.mat'));
                if(max(abs(beta(:) - beta_pre.beta(:))) > 1e-8)
                    error('[Regression from RSFC]: beta differred from original elastic net results')
                end

                feat_test = CBIG_regress_X_from_y_test(feat_test', ...
                    cov_X(test_ind, :), beta);
                feat_train = feat_train';  feat_test = feat_test';
            end

            %% load optimal parameters selected in inner-loops
            params_file = fullfile(LRR_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'params', ...
                ['fold_' num2str(f)], ['selected_parameters_' bhvr_nm{b} '.mat']);
            opt_params = load(params_file)

            %% multilevel block permute y, run elastic net based on permuted y
            for p = 2:Nperm+1
                rng('default')
                rng(1)

                % permute y
                y_perm = y_reg.y_resid(Pset(:,p));
                y_train = y_perm(train_ind);
                y_test = y_perm(test_ind);

                % select features
                if opt_params.curr_threshold ~= 1
                    [feat_train, feat_test] = CBIG_FC_FeatSel( feat_train, feat_test, y_train, ...
                        opt_params.curr_threshold );
                end

                Mdl = fitrlinear(feat_train, y_train, 'ObservationsIn', 'columns', 'Lambda', ...
                    opt_params.curr_lambda, 'Learner', 'leastsquares', 'Regularization','ridge');
                y_p = predict(Mdl, feat_test, 'ObservationsIn', 'columns');

                for k = 1:length(metrics)
                    stats_perm.(metrics{k})(f,p-1) = ...
                        CBIG_compute_prediction_acc_and_loss(y_p, y_test, metrics{k}, y_train);
                end
            end
        end
        save(acc_out, 'stats_perm')
        system(sprintf('ll %s', acc_out))
    end
    
end

function out = reshape3d_to_2d(features)
    % reshapes a #ROI x #ROI x #subjects matrix into
    % #ROI x #subjects by extracting the lower triangle
    % of the correlation
    temp = ones(size(features,1), size(features,2));
    tril_ind = tril(temp, -1);
    features_reshaped = reshape(features, size(features,1)*size(features,2), size(features, 3));
    out = features_reshaped(tril_ind==1, :);
end